Phase-field approach for faceted solidification 



Jean-Marc Debierre^, Alain Karma^, Franck Celestini'^ and Rahma Guerin^ 
^ Laboratoire Materiaux et Microelectronique de Provence (UMR 6137), Universite d'Aix-marseille III, Faculte des Sciences et 

Techniques de Saint- Jerome, Case 151, 13397 Marseille Cedex 20, FRANCE 
^Department of Physics and Center for Interdisciplinary Research on Complex systems, Northeastern University, Boston, MA 

02115, U.S.A. 

Laboratoire de Physique de la Matiere condensee (UMR 6622), Universite de Nice-Sophia-Antipolis, Pare Valrose, 06108 

Nice Cedex 2, FRANCE 

We extend the phase-field approach to model the solidification of faceted materials. Our approach 
consists of using an approximate 7-plot with rounded cusps that can approach arbitrarily closely 
the true 7-plot with sharp cusps that correspond to faceted orientations. The phase-field equations 
are solved in the thin-interface limit with local equilibrium at the solid-liquid interface [A. Karma 
and W.-J. Rappel, Phys. Rev. E53, R3017 (1996)]. The convergence of our approach is first 
demonstrated for equilibrium shapes. The growth of faceted needle crystals in an undercooled melt 
is then studied as a function of undercooling and the cusp amplitude 5 for a 7-plot of the form 
7 = 7o[l + 5{\ sin 9\ + | cos 9\)\. The phase-field results are consistent with the scaling law A ~ 
observed experimentally, where A is the facet length and V is the growth rate. In addition, the 
variation of V and A with 5 is found to be reasonably well predicted by an approximate sharp- 
interface analytical theory that includes capillary effects and assumes circular and parabolic forms 
for the front and trailing rough parts of the needle crystal, respectively. 



I. INTRODUCTION 

Over the last decade, the phase-field approach [1,2] has 
been developed extensively to model the solidification of 
both pure materials and alloys [3] . Most of the work to 
date has focused on the case where the excess free-energy 
of the solid-liquid interface 



7 = 7o/(0) 



(1) 



is a smooth function of the angle 9 between the direction 
normal to the interface and some fixed crystalline axis. 
In particular, the simple form 



/(6») = H-ecos46'- 



(2) 



appropriate for a weakly anisotropic material with an un- 
derlying cubic symmetry has been widely used in studies 
of dendritic solidification [4-9]. For a smooth 7-plot, the 
value of the diffusion field (dimensionless undercooling or 
supersaturation) at the interface is given by the standard 
Gibbs-Thomson condition 



-do 



(3) 



where do is a microscopic capillary length (thermal or 
chemical) that is proportional to 70 and n is the inter- 
face curvature. Dendritic growth in pure materials [6-8] 
and alloys [9] has been modeled quantitatively with a 
thin-interface limit [6] of the phase-field model that yields 
the boundary condition (3) and overcomes the stringent 



computational constraint associated with a finite inter- 
face thickness. These simulation studies modeled the 
common situation of a weak anisotropy where the inter- 
face stiffness 7 -I- cP^/dO"^ , and hence the square bracket 
in Eq. (3), is everywhere positive. More recently, the 
phase-field approach has also been successfully extended 
to model strongly anisotropic equilibrium crystal shapes 
where the stiffness becomes negative for some range of 
9 [10]. In this case, which occurs for example when 
e > 1/15 for the form (2) of f{9), the slope of the in- 
terface has discontinuities that must be present to elim- 
inate thermodynamically unstable, and hence forbidden, 
orientations for which the stiffness is negative. 

In this paper, we extend the phase-field approach to 
model a wide class of materials that form facets for a 
discrete set of orientations. The interfacial energy is gen- 
erally non- analytic for orientations close to a facet. This 
non-analyticity is reflected in the presence of cusps in the 
7-plot that are of the form 



f{9)^\ + 5\9-9,\ + 



for 161 - I 



« 1, 



(4) 



where 9 — 9c. is the orientation of a given facet. In the 
simplest model where vicinal surfaces are assumed to con- 
sist of straight ledges of height a separated by terraces of 
width I « a/\9 — 9c\, the cusp amplitude 5 — Jl/{o,1o) 
where 7^ is the ledge energy per unit length. 

To a good first approximation, the tendency of a ma- 
terial to facet is described by Jackson's well-known a- 
factor [11] that is the product of a dimensionless crys- 
tallographic factor that depends on the crystal structure 
and the orientation of the interface and the ratio of the 



1 



latent heat per mole to the rare gas constant. Jackson's 
theory predicts that materials with a < 2 (such as metals 
with low entropy of melting) will have rough interfaces, 
while those with a > 2 (such as non-metals and com- 
pounds) will form facets, consistent with experimental 
observations across a wide range of materials [12]. 

The growth of faceted dendrites has been studied both 
experimentally [13,14] and theoretically in the context 
of a sharp-interface approach [15 -17]. Experimentally, 
facets are seen to appear near the tip of the needle crys- 
tal that nevertheless retains a parabolic shape on the 
larger scale of the whole crystal that includes the trailing 
rough parts. From a theoretical standpoint, a main diffi- 
culty to model faceted growth is that the standard Gibbs- 
Thomson relation can only be applied to rough parts of 
the interface. On a facet, this relation takes the form 
of an integral condition that is obtained by making the 
substitution k = d6/ds, where s denotes the arclength 
along the interface, and integrating both sides of Eq. (3) 
from one extremity (s_) of the facet to the other (s+). 
This integration yields the condition 



/' 

J S- 



uds 




d0 



-2 do 5, (5) 



where df/dd\± denotes the limiting values of df/d6 on 
each side of the cusp {0 — 9c 0"^). 

For an equilibrium crystal with a constant undercool- 
ing or supcrsaturation (u < 0) along the interface, the 
above condition implies that the length A of the facet is 
simply proportional to the amplitude of the cusp 



—u 



(6) 



For growth outside of equilibrium, the determination of 
the facet length and the shape of the rough parts requires 
in general a self-consistent solution of the free-boundary 
problem defined by the appropriate diffusion equations 
for u in solid and liquid, the standard Stefan condition 
of heat or mass conservation on the interface, and the 
local equilibrium conditions (3) and (5) on rough parts 
and facets, respectively. We defer a discussion of kinetic 
effects to the conclusion section of this paper. 

An analytical solution to this problem was obtained for 
steady-state growth by Adda-Bedia and Hakim [16] in the 
small Peclet number limit neglecting capillary effects on 
the rough parts. They concluded from this analysis that 
it is not possible to require both tangential matching of 
the rough and faceted parts of the interface and the equi- 
librium condition [u « 0) on the front and trailing rough 
parts. Furthermore, they proposed an approximate so- 
lution to the full problem that includes capillary effects 
on the rough parts in the large 6 limit. In this solution, 
the facets are assumed to match tangentially to a small 
quasi-circiilar tip that is significantly more undercooled 
than the rough trailing parts. A numerical solution of 
the steady-state growth problem was later obtained for 



arbitrary 5 by Adda-Bedia and Ben Amar [17] using a 
boundary integral method. This calculation gave results 
that are consistent with the the approximate solution of 
Adda-Bedia and Hakim for large 5. 

At present, a phase- field approach for faceted growth 
would be highly desirable to study the stability of needle 
crystal solutions as well as to explore the full dynami- 
cal range of morphological evolution when the interface 
shape is non-stationary. Treating separately the rough 
and faceted parts of the interface appears to be difficult 
within a phase-field formulation. To circumvent this dif- 
ficulty, we follow a procedure that is similar in spirit, 
but different in details, to the one introduced by Adda- 
Bedia and Ben Amar [17] in a sharp-interface context. 
The basic idea is to use a regularized form of the 7-plot 
that can approximate arbitrarily closely the 7-plot with 
sharp cusps. The discontinuity of df/dd at a cusp can 
be viewed mathematically as the stiffness behaving as a 
sharply peaked Dirac delta function at ^ = ^c- There- 
fore, different regularizations of this discontinuity can be 
viewed as different regularizations of this delta function 
peak that must satisiy the cusp condition 



lim 



i 



dO 



^ de^ 



26, 



(7) 



imposed by the form (4) of / near a cusp. Adda-Bedia 
and Ben Amar have introduced regularizations in which 
all derivatives of the stiffness and hence / are continu- 
ous. In a phase- field context, we have found more con- 
venient to use a regularization where only / and its first 
derivative are continuous and the dimensionless stiffness 
/ -|- d^f/dO"^ is a step function regularization of the delta 
function, f + d'^f/dO'^ « (5/6*0 for 9c-0q <9< 9c + 0o- In 
practice, computations for small ^0 turn out to be more 
costly because the larger stiffness (« S/60) requires a 
smaller time step. As we shall see, however, an accurate 
extrapolation to the ^0 — * limit is nevertheless possible. 

The phase-field model and the details of this smoothing 
procedure are described in the next section. The numer- 
ical implementation of the equations is then presented in 
section HI. The convergence of our numerical approach 
is demonstrated in section IV by a detailed comparison 
of diffuse-interface and sharp-interface equilibrium crys- 
tal shapes. Results concerning the steady-state growth 
of faceted needle crystals in a pure undercooled melt are 
then presented in section V. These results are then com- 
pared in section VI to the predictions of an approximate 
sharp-interface theory that assumes simple forms for the 
front and trailing rough parts as an extension of the so- 
lution proposed in Ref. [16]. Concluding remarks and 
future prospects are then discussed in section VI. 



II. PHASE-FIELD MODEL 

We consider the growth of a crystal from a pure under- 
cooled melt with the thermal diffusivity being the same 
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in the solid and liquid. Wc follow the phase-field ap- 
proach developed by Karma and Rappcl [6] that has 
proven successful to model cfBcicntly non-faceted den- 
dritic growth. We briefly recall the basic features of the 
phase-field model and of the thin-interface limit used to 
relate the phase-field and sharp-interface models. 

The dimensionless temperature field is defined as u = 
(T — Tm)I{L/cp), where T is the temperature, Tm the 
melting temperature, L the latent heat per unit volume, 
and Cp the specific heat at constant pressure per unit 
volume. The phase field V) is taken to be equal to +1 
(— 1) in the solid (liquid) phase, and varies continuously 
across the diffuse interface. In two dimensions, the time 
evolution of both fields is governed by the equations 



dtu = DV^u + i^tV, 



(8) 



and 



[V'-An(i-V')](i-V'') 

+v- [w{6fvij] -d^[w{e)w'{e)dy7p] 

+dy[w{e)w'{e)d,ij], (9) 



where D is the thermal diffusivity, 9 is the angle between 
the normal to the interface and the x axis, A is a dimen- 
sionless coupling constant between the temperature and 
phase fields, and 



wie) = Wofie), 



(10) 



is the diffuse interface thickness. In the thin-interface 
limit where Wq is much smaller than the mesoscale of 
the growing crystal [6], the above phase- field equations 
reduce to the standard sharp- interface model of diffusion- 
limited growth with the velocity-dependent form of the 
Gibbs-Thomson condition 



-do 



where 



and 



m + 



do 



d'm 

de^ 



aiWo 



m = 



ai t{9) 



A w{e) 



1 — a2A 



Dt{6) 



(11) 



(12) 



(13) 



In addition, ai and a2 are constants that depend gen- 
erally on the choice of the double-well potential and 
other functions in the phase-field model. For the present 
choices that are the same as in Refs. [6], a\ = 0.8839 . . . 
and 02 = 0.6267 . . .. 

To model faceted growth, we use the simplest form of 
the 7-plot 



f{0) = l + d{\sme\ + \cose\). 



(14) 



This form is directly adapted from the broken-bond 
model which describes satisfactorily faceted solid-gas in- 
terfaces [18]. The same form was used in previous sharp- 
interface calculations [16,17]. 

Furthermore, in order to model the limit with local 
equilibrium at the interface, i.e vanishing interface kinet- 
ics (}{9) = 0, we impose [6] 



r{9) = Tof{9f 



and 



A = 



1 Dto 
a2 Wi ' 



(15) 



(16) 



which makes P{9) vanish for all values of 9. 

For the above choice of 7-plot, the dimensionless stiff- 
ness defined by 



S{9) ^ f{9) 



d9^ 



(17) 



is constant except at the cusps, 9 = 9c,n = mr/i (with 
n an integer), where it diverges. The discontinuities in 
the first derivative of the interfacial energy need to be 
regularized, since f'{9) is not defined for 9 = 9c,n- A 
simple strategy to circumvent this problem is to smooth 
out the cusps by replacing f{9) with a smooth function 
like a sine in a small range of 9 values around the cusps 
(Fig. 1). For the sake of clarity, we restrict our discussion 
to the first quadrant, 9 G [0, §], since the problem has 
a natural four-fold symmetry. The smoothed anisotropy 
function fs{9) reads then: 

1 -I- (5 (sin 6* -I- cos ^) 

f^{e) = { B-Acos9 
B - Asm9 



if 

if 
if 



^0 < 9 < ^ — 9o, 

9 < On. 



(18) 



and its derivative 

(5(cos^ — sm9) 
m ={ Asm9 
-Acos9 



if 6*0 < 6* < f 
ff 6* < 6*0, 
if 9>l-9o. 



(19) 

The two constants A and B are obtained by expressing 
the continuity of fs and f'^ at 9 = 9o, 



and 



A = 5(cot(9o - 1), 



B = l + S/sm9o. 



(20) 



(21) 



Now /s(0) = /s(7r/2) = and W{9) is no more singular 
in the cusp directions. Note, however, that f'J is not 
continuous at ^ = ^o- Consequently, the dimensionless 
stiffness is a step function of 9 (see Fig. 1), 



3(9) = 



1 



1 + 5/ sin 6*0 



if 
if 



9o<9<^ 
0<0o or 



-9o, 
9>l 



60. 
(22) 
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III. NUMERICAL IMPLEMENTATION 

We now briefly describe how Eqs. (8, 9) are discretized 
in our code. The interface is represented by the ■0 = 
contour, so that V0 is by construction cohnear to the 
unit vector along the normal to the interface, n. If 
I VV'I 7^ 0, the two components of this vector are given by 



= cos^ = -a^^/|VV'|, 



and 



sm f 



(23) 



(24) 



and fs and are computed according to Eqs (18, 19). 
Conversely, when [VV'I = 0, one sets 



and 



fs{0) 



(25) 



(26) 



Numerical simulations of Eqs. (8, 9) are performed by im- 
plementing a finite-difference scheme on a regular square 
mesh with mesh size h = OAWq. The domain considered 
is the quadrant < .t < Nh and < y < Nh, with an 
integer. A standard first-order in time Euler scheme with 
a time step At and a second-order in space discretization 
of the spatial derivatives are used. 

Let Ui,j and respectively denote the discretized re- 
duced temperature and phase field at point {x = ih, y = 
jh). In a first step, one computes V^u, V^tp, d^fp, dytp, 
W{9)W'{e) and W{e) for each point of the extended do- 
main, i and j € [— 1, A'' -(- 1], and the results are stored 
in six intermediate arrays. Reflecting conditions are im- 
posed on u and tp at all the domain boundaries {i or 
j = —2,—l,N+l,N + 2), and using centered differenc- 
ing approximations, we have 

(27) 



0. 



01, 



i-1 



dxfpij = [tpi+ij -i>i-i,j]/{^h), 
dyipij = [tpij+i - V'i,j-i]/(2/i)- 



(28) 



(29) 



(30) 



The two remaining arrays, (WW')ij and Wij, are given 
by Eq. (10) into which / is replaced with /g. 

In the second step, one solves Eqs. (8) and (9) on the 
inner points of the domain, i and j G [0, A^]. To do so, 
the last three terms in Eq. (9) are written in a slightly 
different form, 



W{efV^^ + 2W{9)VW(6)Vtlj 

-{dy^)d,[w{9)w'{9)] + {d,^p)dy[wi9)w'i 



(31) 



so that each new term can be computed with the help 
of the intermediate arrays. For instance the last term is 
discretized as 

(5.Vi,i) [( W)i,,-+i - {WW%,.^] /{2h). (32) 



IV. EQUILIBRIUM SHAPES 

A. Comparison between analytical and phsise-fleld 
results 

We first check our code by computing the equilibrium 
shape of the crystal. In thermal equilibrium, u = —A 
everywhere, and the evolution equation for the phase field 
reduces to 

T{9)dt^ = + AA(1 - V'^)] (1 - + w{efv^^ 

+2W{9)VW{9)V^j: - {dy'4))d^ \W{9)W'{9)\ 
+{d,i,)dy[W{9)W'{9)\. (33) 

Initially, the solid is given a circular (or a square) shape 

and the undercooling is set to A = Aq. Time integration 
of Eq. (33) is performed and the interface velocity along 
the X axis, 14 is computed at regular time intervals. If 
\4 > 0, A is decreased and it is increased if 14 < 0. The 
increment JA is divided by a constant n > 1 each time 
Vx changes sign. This procedure is rc^peated as long as 
5 A is larger than some prescribed value, dAmin- This 
scheme is known to converge to an equilibrium state [6] . 
The analytical equilibrium shape is given by 



I' 

I y 



X = x{9)/Ro 
y = y{9)/Ro 



- fs{9) cos 9- f'M sin 9, 
/,(e)sin^ + /^(e)cos^, 



(34) 



where Rq = do /A [19] . Let us consider one eighth of this 
interface corresponding to 6* G [7r/4,7r/2] (Fig. 2). In the 
cusp region, 9 G [n/2 — 9o,Tr/2], one has 



x{9) =Bcos9, 
y{9) = {B sin 9 -A), 



(35) 



so that the interface is the circle of center (0, —A) and 
radius B. Since the smoothed cusps are very narrow in 
practice (^o 1) this circle arc is very flat; it tends to 
the straight horizontal facet as — * 0. The right end of 
this interface portion lies at 



.XV = x{7r/2 - 9o) = Bsin9o, 
yr = y(7r/2 - 6»o) = B cos 6*0 - A, 



where the local slope is 



■ tan^o- 



(36) 



(37) 
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On the other hand, for 6 G [7r/4, 7r/2 — 6*0] the crystal 
is bounded by the circle centered at point Pc{5,6) (i.e. 
X = y = 5) of rescaled radius unity, for which 



R={i^- 5^)V2 = 1 - \/2sin6lo 



(38) 



is a good first-order approximation when ^ 1 (Fig- 2). 

Note that, as 9q 0, one recovers the equilibrium 
shape for the sharp cusp, with a horizontal facet of length 



(39) 



matching tangentially a circle; of rescaled radius unity. 

As shown in Fig. 3, the whole crystal shape is well 
reproduced by the phase-field equilibrium code. A close 
examination of the nunieric:al data points reveals that 
the imposed anisotropy is underestimated by about 0.2 
percent in the numerics, independently of ^o- 



B. Numerical estimates for the facet ends 

The knowledge of the analytical equilibrium shape 
guides us to define a numerical procedure to extract es- 
timates of the facet length A and corner radius R from 
discrete interfaces obtained with the phase-field code. On 
the square mesh, the interface consists in a list of points 
Pi{xi,yi)- According to Eq. (37), the two facet ends are 
located at the points Pr and Pi where the absolute value 
of the local slope rapidly increases beyond tan^?o (Fig. 4). 
The method used to compute the two end points is thus 
to compare the derivative of the equilibrium curve dxU 
with tan^o- Using centered differences. 



dxVi = 



Ui+i - Ui-i 



(40) 



gives an error of about 5 percent on both measures. On 
the other hand, using the one-sided approximations, 



dxVil 

to compute Pr and 



Vi - Vi-i 



Vi+i - Vi 

Xi-\-'l Xi 



(41) 



(42) 



to compute P/, estimates of A and R are only 0.5 percent 
off the exact values (see Table I ). This last estimation 
procedure is thus more precise and it is readily extended 
to the case of growth shapes in what follows. 



T{e)dti> = [v- + AA(i - V/)] (1 - V') 

+ Wq [Cxxdxxi^ + CxydxyIP + Cyydyytp] , (43) 



with 



1 + (5(3|cos6l| - Icos^e*! + \sm^e\)+S^, 

- AB\ cos 0\{2 + sm^0), (44) 
B'^ - AB\sm^0\, 



2S{as cos^ 9 + ac sin^ 6) + 2S^GsOc 



Cxv = { -2AB(jc 



and 



yy 



-2ABas cos^ 9, 



1 + 5(3| sin 6*1 - I sin^ 6*1 + \ cos^ 6*1) + 5^ 

B"^ - AB\cos^e\, 

A^ + B^- AB\ sin 6*1 (2 + cos^ 6), 



(45) 



(46) 



hr 9o < 9 < ^ - 9f), 9 < 9o, and 9 > f - ^o, respectively. 
The symbols as and Uc represent the signs (±1) of sin 9 
and cos^. 

The time step At must be siifficicntly small to en- 
sure convergence of the finite difference scheme, the most 
stringent constraint on At arising for 6 = Oq: 



At < 



Wi 2 



Cxx{6o) + Cyy{6o) 



As ^0 ^ 1; this approximates to 



(47) 



(48) 



It is thus always possible to reduce 9q provided that At 
is decreased in proportion. Although At depends on the 
cusp amplitude i5, the time step can still be kept constant 
while maintaining numerical stability by imposing 



9^ = K 



1 + 6' 



(49) 



with K some constant. In the numerical simulations, we 
tipically take K = tt/IOO, so that At/ro < 0.0025.... 
In practice, we rather use the discretization scheme de- 
scribed in section III than the fully developped version 
given in Eqs. (43-46). The former proves to be more pre- 
cise and more stable than the latter, so that a larger time 
step, At/ro = 0.008, can be used. 



V. FACETED NEEDLE GROWTH 



C. Sharp cusp limit 

One may still wonder if the sharp cusp limit (^o — * 0) 
is reacheablc within our numerical approach. To answer 
this question, it is better to develop Eq. (33) in terms of 
the second partial derivatives of One gets then 



We now turn to non-equilibrium growth shapes. The 
time evolution of a needle dendrite is illustated in 
Fig. 5 for a cusp amplitude ^ = 1, diffusion coefficient 
Dt^/Wq = 4 and imdercooling A = 0.55. The size of the 
simulation box is 6001^0 x 600 Wq. For reasons of symme- 
try, it is sufficient to grow one half of the dendrite (here 
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y > x). We checked that identical patterns are obtained 
when the whole quadrant is used. The initial conditions 
for the phase field are ip = 1 ii both < a; < 20Wo and 
< y < 20 Wo and tp = —1 otherwise. This square shape 
is imposed to the solid germ in order to get a dendrite 
with a single tip along the x = y direction. Conversely, 
rounded initial shapes, like a circle, always result into 
tip splitting along the growth direction, leading to more 
complicated patterns which will be discussed later. The 
temperature field u is initially set to zero inside the solid 
germ and to —A outside. Reflecting conditions are im- 
posed at all the domain boundaries 

We first focus on the influence of the anisotropy co- 
efficient, (5, on the dynamics of the faceted needle den- 
drite. This coefficient is varied from 0.25 to 1.60 in the 
simulations. After a transient, the dendrite reaches a 
stationary state which is independent of the size of the 
initial germ. The tip then moves with a constant veloc- 
ity V and a rather large portion behind the tip adopts 
a stable shape. A good test of the thin-interface limit is 
to decrease independently the rescaled interface thickness 
Wo/rfo (or equivalently the diffusion coefficient Dtq/Wq) 
and the grid spacing h/Wo [6]. Convergence is achieved 
when the scaled tip velocity V = Vdo/D is independent 
of both pameters. Table II shows that the results are rea- 
sonably well converged for the values Dtq/Wq = 4 and 
h/Wo = 0.4 used in our simulations. 



A. Steady-state shape 

The steady-state interface shape can be divided into 
three distinct parts: a nearly circular tip of radius R, 
a smooth "facet" of length A, and a trailing rough tail 
(Fig. 6). Estimates of R= [xt— Xr)V^ and A = Xr — xi 
are obtained by using the numerical interpolation method 
described in section IV-B. To smooth out temporal fluc- 
tuations, time averages of both quantities are performed 
over the whole stationary regime. Obviously, the exact 
shape of the dendrite tip should depend on the smoothing 
angle 6o. When Oq is sufficiently small this dependence 
is found to be linear, which allows us to unambiguously 
extrapolate R{9o) and A{9o) to 6*0 = (Fig. 7). Up to an 
overall scale factor, the dendrite tip shape is very sim- 
ilar to the equilibrium shape, i.e. a quarter of a circle 
matching tangentially to two side facets. 

The variations of A(0) and -R(O) with cusp amplitude 
are dispayed in Fig. 8. The results are well fitted by the 
two simple relations. 



and 



A = Ao 



R = Rq 



1 



S 



Sr 
S 



(50) 



(51) 



The asymptotic shape of the trailing rough part away 
from the facet is well-fitted by a parabola that satisfies 
the Ivantsov relation [20] 



with the Peclet number 



(52) 



(53) 



Note that the reason why the asymptotic tail is a 
parabola whose tip radius is predicted exactly by the 
Ivantsov relation is the same as for non-faceted growth 
[21,22]. Namely, in a frame that is stationary with re- 
spect to the melt, the trailing rough part grows as a 
planar interface whose position x{t) perpendicular to 
the growth axis ^ t^^"^ and the Peclet number p = 
x{t){dx{t)/dt)/2D obeys the same relation as Eq. (52), 
as shown independently of Ivantsov by Zener [23]. 

We checked in a few simulations that pV is well 
predicted by the Ivantsov relation, as illustrated in 
Fig. 9 where we compare the; phase-field interface to 
the Ivantsov parabola with the same tip velocity, i.e., 
the parabola with p = 2Dp/V, where p is the Ivantsov 
Peclet number and V is the steady-state tip velocity in 
the phase-field simulation. Given this agreement, we use 
directly the formula p = 2Dp/V to calculate p in all re- 
sults reported thereafter. It should be emphasized the 
tip radius p of this imaginary parabola that matches ex- 
actly the asymptotic needle crystal shape and the true 
tip radius R are quite different as shown in the plot of 
R/ p versus 3 in Fig. 10. In particular, R/ p is seen to de- 
crease sharply with increasing 5, while A/p approaches a 
constant. These plots reflect the fact that side facets be- 
come elongated and extend closer to the tip that becomes 
more pointed as 5 increases. 



B. Steady-state operating state 

As in previous sharp-interface calculations [16,17], we 
define the dimentionless tip selection parameter 



C 



Ddo ' 



(54) 



where V is the tip velocity and p is the tip radius of 
the parabola that exactly matches the asymptotic tail of 
the faceted needle crystal. Note that C = 8/a where 
a = 2Ddo/p'^V is the selection parameter that is also 
often used in dendrite growth studies. 

Eliminating p in favor of the Peclet number using the 
relation p = 2Dp/V , we can rewrite C in the form 



C ■ 



IQD 

Vdo 



P ■ 



(55) 



A plot of C versus 6 with V corresponding to the steady- 
state tip velocity in the phase-field simulations and p 
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computed from the Ivantsov relation for the value of A 
corresponding to the simulations is shown in Fig. 11. 
Note that C{5) has a minimum that corresponds to a 
maximum of tip velocity as a function of 5. 

In experiments, the growth rate of a single material is 
usually studied as a function of undercooling Additional 
simulations are thus performed for other values of A. All 
the other parameters are kept constant, except for the 
total time which is increased to TOOOtq for A = 0.45, 
3000ro for A = 0.50, and reduced to 1500to for A = 
0.60. We do not extrapolate V to 6q ^ here because 
the computing time becomes prohibitively large for A < 
0.55. However, from the previous results for A — 0.55, 
we do not anticipate deviations in the tip velocity by 
more than a few percent (see Fig. 7). The plot of tip 
velocity versus facet length shown in Fig. 12 shows that 
our numerical results are consistent with the scaling law 
A ~ V~0-5 which was found experimentally for NH4Br 
needle crystals [14]. 

C. Transients 

When growth is started from a circular germ, a dou- 
ble needle is systematically observed (Fig. 13), indepen- 
dently of the germ shape, for 5 < 0.20. Initially, the two 
tips move away from each other. Our simulations show 
that these needles ultimately grow practically parallel to 
each other at very long times, wh(ui they arc separated 
from each other by several diffusion lengths (Fig. 13). 

At long times, the tip of the double needle becomes 
comparable in shape to that of the single needle dendrite 
(Fig. 14). Of course, the former is thinner on one side 
because the latent heat accumulates in the channel be- 
tween the two tips. This similarity in shape suggests that 
the double needle pattern is governed by the same oper- 
ating state as the needle dendrite. Verifying this point 
with a sufficient accuracy is not easy because of a still 
longer transient for the double needle dendrite. For this 
reason, we compare the patterns obtained for a larger un- 
dercooling, A — 0.60, for which the simulations are still 
quantitative but the relaxation time much shorter than 
with A = 0.55. After a time t = 2000ro, the tip velocity 
of the double needle pattern has not yet fully converged 
but it tends slowly to a limit near 0.397VFo/to (Fig. 15). 
This is precisely the value found for the single tip den- 
drite, which confirms that the two asymptotic states are 
identical. Thus, the double needle dendrite is merely just 
a tip-splitting evolution of the single tip dendrite. This 
splitting is often preferred by the system because it allows 
for a more efficient occupation of the available room. 

Let us finally remark that the double needle found here 
is not a faceted version of the generic non-faceted growth 
structure that has been termed doublon [24]. Doublons 
also are symmetry-broken growth shapes with a double 
tip. However, the two tips remain much closer than 
one diffusion length and are thus highly inter-dependent. 



which is not true for the double faceted needle dendrite 
observed in the present simulations. 

VI. COMPARISON WITH SHARP-INTERFACE 
ANALYTICAL THEORY 

A precise quantitative comparison of the phase-field 
results of the last section with previous sharp-interface 
calculations [16,17] cannot be carried out because these 
calculations considered the one-sided model with zero dif- 
fusivity in the solid-phase while we simulated here a sym- 
metric model with equal diffusivities in solid and liquid. 
In addition, our simulations are for finite Peclet num- 
ber while the sharp-interface calculations of Refs. [16,17] 
are restricted to the small Peclet number limit (p — ^ 0). 
Nonetheless, the phase- field results reproduce several key 
features of faceted needle growth that were predicted in 
the sharp-interface studies: 

1. The ratio A/p of the facet length to the tip radius 
of the parabola that matches the asymptotic rough 
tail increases with the cusp amplitude S and satu- 
rates to a value of order unity for large S. 

2. The tip selection constant C has a minimum for 
a value of 5 of order unity. The value of C that 
corresponds to this minimum is large 100). 

3. For large S, the needle shape consists of two long 
facets joined by a nearly circular tip that match on 
to trailing rough parts that become parabolic far 
from the tip. This is precisely the shape proposed 
by Adda Bedia and Hakim [16]. 

In the remainder of this section, we develop a simple 
analytical theory of faceted needle growth that is based 
on the approximate shape proposed by Adda Bedia and 
Hakim [16] for large S, where facets are joined by a small 
circular tip of radius R. We take this approximation fur- 
ther by assuming that (i) it remains valid for 6 of order 
unity, and (ii) the trailing rough parts remain parabolic 
all the way to the points at which they match tangentially 
the facets, as opposed to being parabolic only asymptot- 
ically far from the tip [16]. These simplifications allow 
us to obtain a simple physical picture of faceted needle 
growth as well as explicit analytical predictions of the 
shape parameters R, A and p and the velocity V for ar- 
bitrary Peclet number, without necessitating the numer- 
ical solution of an integral equation as in Refs. [16,17]. 
As \vc shall sc!e below, thcisc; prcidictions agree reason- 
ably well with our phase-field results despite our simple 
parametrization of the steady-state shape. 

The complete determination of the needle shape {R, 
A and p) and the growth velocity V for fixed A and 
6 requires four independent relations. The first is the 
Ivantsov relation (52) that fixes the product pV. A sec- 
ond relation between R, A and p is simply obtained by 
imposing that the circular tip and parabolic tails match 
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tangentially the front and trailing ends of the facets, re- 
spectively, which yields at the relation 



p = (A + R)/V2. 



(56) 



We can check this relation by comparing the tip ve- 
locity in the phase-field simulations with the velocity 
V = 2Dp/ p, where p is the Peclet number predicted by 
the Ivantsov relation and p is computed using Eq. (56) 
together with values of A and R in the simulations. 
Fig. 16 shows that these two velocities are in reasonably 
good quantitative agreement. The systematic deviation 
is most likely due to the fact that the crystal tail is not 
exactly parabolic very close to the facet. 

The two other relations needed to complete our theory 
are obtained by imposing the Gibbs-Thomson condition 
(3) at the tip of the needle crystal and by using the in- 
tegral form of this condition on the facet (Eq. 5). For 
this, we need an expression for the undercooling along 
the interface that is obtained straightforwardly using the 
known boundary integral relation [25,26] 

/•+°° dx' 

u(a;) = — A + p / — eici){—p[z{x) — z{x')]) 

J~oo I" 

xKq {p^{x-x'y + {z{x)-z{x')Y'^ , (57) 

where Kq is the zeroth order modified Bessel function, 
p is given by the Ivantsov relation (52), and z{x) is the 
interface shape with length in unit of p. For a circular 
tip and parabolic tails matching tangentially side facets 
that make a 45° angle with respect to the growth axis, 
the X coordinates of the front and trailing ends of the 
facets are x = ±R/V2p and x = ±1, and 



z = 



R^_ 

-\x\ + 

_x^ 
2 ' 



1/2 



R 



X 

1 

2' x/2p 

\x\ > 1. 



V2R 1 

< \x\ < 1, 



Q<\x\< 



R 



(58) 
(59) 
(60) 



These equations describe rather accurately the phase- 
field needle dendrites (Fig. 17). Since the isothermal 
Ivantsov parabola is an exact solution of the steady-state 
growth problem without capillarity, the integral relation 
(57) is exactly satisfied for z = —x^/2 and u = 0. In 
contrast, as illustrated in Fig. 9, the tip of the faceted 
needle crystal protrudes ahead of the tip of this parabola 
and hence is undercooled by a finite amount 'u(O). The 
Gibbs-Thomson condition (3) imposes a relation between 
this tip undercooling and the tip radius that is simply 



u{0) = -do/R. 



(61) 



In turn, Eq. (5) provides us with a relation between the 
average undercooling on the facet and the cusp amplitude 
S that takes the form here 



V2p 



dxu(x) = —2— i5, 
P 



(62) 



where x is in units of p as above. The four relations de- 
fined by Eqs. (52, 56, 61) and (62), together with the 
expression for u on the interface defined by Eqs. (57-60), 
completely determine R, A, p and V. Since the product 
pV is exactly predicted by the Ivantsov relation, we only 
need to compute three independent dimensionless combi- 
nations of the four above quantities to compare the pre- 
dictions of the above theory with the phase-field results. 
The most meaningful dimensionless combinations are the 
selection parameter C and the two ratios A/p and R/p. 
To compute those as a function of 6, it is convenient to 
vary R/p and compute C, 5, and A/p using the following 
relations that are simple to deduce: 



1 _ 
C ' 

s = 



V2 

16pS I R 

V2p 



dxu{x), 



^/2u{0)R , 



dxu{x), 



s/2p 

R 
72^ 



(63) 
(64) 
(65) 



The above predictions are compared with the phase-field 
results in Figs. 10 and 11. The quantitative agreement 
between the two is reasonably good but not exact because 
the steady-state shape that we have assumed here difl:ers 
slightly from the true shape (Fig. 17) that is only exactly 
parabolic asymptotically far from the facet. 

To conclude, we note that the expression for u{x) can 
be simplified in the small Peclet number limit where 
Eq. (57) reduces to [27] 



u{x) 



dx^ 
7-00 27r 



In 



{x - x')^ + {z{x) - z{x'yf 



{x-x')^ + {x'^/2-x^/2)^ 



(66) 



Therefore, both w(0) and the integral of u{x) along the 
facet are simply proportional to the Peclet number in 
this limit and hence C, R/p, and A/p are only functions 
of d. Furthermore, it is straightforward to deduce that 
A/p ^ \/2, R/p 0, and C ~ (5 in the large S limit. 
Hence this theory predicts the scaling law A ^ 
that is a simple consequence of the fact that A/p and 
C ~ p'^V are constants in the small Peclet number limit. 



VII. CONCLUDING REMARKS 

In conclusion, we have shown that the phase-field ap- 
proach can be successfully extended to model the solid- 
ification of faceted materials. Our approach, which con- 
sists of rounding the cusps in the 7-plot, converges well 
in the limit of sharp cusps for both equilibrium and non- 
equilibrium growth shapes with facets. Even though we 
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have considered a simple form of the 7-plot, the method 
should be applicable to more complex 7-plots where the 
interface stiffness varies on the rough parts. 

In addition, we have developed an approximate ana- 
lytical theory of faceted needle growth that includes cap- 
illarity and assumes circular and parabolic forms for the 
front and trailing rough parts of steady-state needle crys- 
tals, respectively. This theory yields explicit predictions 
of tip velocity and facet length that arc in good overall 
quantitative agreement with the phase-field results; this 
agreement is largely due to the fact that the needle shape 
assumed in the theory is a very good approximation to 
the actual growth shape observed in the phase-field sim- 
ulations. Furthermore, this theory predicts that the scal- 
ing law A ~ V~^/^ observed both experimentally [14] and 
in the present phase-field simulations should hold in the 
small Pcclet number limit, consistent with the results of 
previous theoretical studies [15-17]. 

In the present study, we have established a quantitative 
comparison between phase-field and sharp-interface re- 
sults under the assumption that the variation of interface 
undercooling is dominated by capillarity. Using scaling 
arguments, Ben Amar and Pomcau [15] have concluded 
that, under the assumption that the growth of facets is 
dominated by a Franck-Read screw discolation mecha- 
nism, the kinetic undercooling of the interface should be 
negligible over an intermediate range of velocity where 
the A ~ scaling law has been observed [14]. While 

growth conditions where kinetic effects are negligible may 
exist, we suspect that such effects will generally be impor- 
tant in the presence of facets, as widely believed. There- 
fore, a quantitative incorporation of facet kinetics in a 
phase-field model remains an important task for the fu- 
ture. The phase-field model as formulated thus far re- 
produces a linear relationship between the planar inter- 
face velocity and interface undercooling appropriate for 
a rough interface. On facets, however, standard growth 
mechanisms such as screw dislocations or ledge nucle- 
ation lead to a non-linear relationship between velocity 
and undercooling. These relationships could potentially 
be incorporated in the phase-field model by letting the 
kinetic relaxation time t depend on temperature or su- 
persaturation, in addition to orientation. 

Another interesting future prospect is to model the di- 
rectional solidification of alloys with faceted interfaces 
[28-30] by combining the present methodology to han- 
dle cuspcd 7-plots with a recent thin-interface phase-field 
formulation of alloy solidification [9] . 

ACKNOWLEDGMENTS 

A.K. wishes to thank the hospitality of the Labora- 
toire Materiaux et Microelectronique de Provence, Uni- 
versite d'Aix-marseille III where this work was initiated 
as well as the support of U.S. DOE Grant No. DE-FG02- 
92ER45471. We also thank Vincent Hakim and Klaus 



Kassner for valuable discussions. This research benefited 
from computer time allocation at the Advanced Scientific 
Computation Center at Northeastern University. 



[1] J. S. Langer, in Directions in Condensed Matter (World 

Scientific, Singapopre, 1986), p. 164. 
[2] J. B. Collins and H. Levine, Phys. Rev. B 31, 6119 

(1985). 

[3] For recent reviews, see W. J. Boettinger, ,J. A. Warren, 
C. Beckermann, and A. Karma, Ann. Rev. Mater. Res. 
32, 163-194 (2002); A. Karma, "Phase Field Methods", 
in Encyclopedia of Materials: Science and Technology, 
edited by K. H. J. Buschow, R. W. Cahn, M. C. Flem- 
ings, B. Ilschner, E. J. Kramer, S. Mahajan, Volume 7, 
Elsevier, Oxford, pp. 6873-86 (2001). 

[4] A. A. Wheeler, B. T. Murray, and R. Schaefer, Physica 
D 66, 243 (1993). 

[5] S-L. Wang and R. F. Sekerka, Phys. Rev. E53, 3760 
(1996). 

[6] A. Karma and W. J. Rappel, Phys. Rev. E 51, R3017 

(1996); Phys. Rev. E 57, 4323 (1998). 
[7] N. Provatas, N. Goldenfeld, and J. Dantzig, Phys. Rev. 

Lett. 80, 3308 (1998); ,L Comp. Phys. 148, 265 (1999). 
[8] M. Plapp and A. Karma, Phys. Rev. Lett. 84, 1740 

(2000); J. Comp. Phys. 165, 592 (2000). 
[9] A. Karma, Phys. Rev. Lett. 11, 115701 (2001). 
[10] J. J. Eggleston, G. B. McFadden, and P. W. Voorhees, 

Physica D 150, 91 (2001). 
[11] K. A. Jackson, in Liquid Metals and Solidification (ASM, 

Cleveland, Ohio) p. 174. 
[12] D. P. Woodruff, The Solid-Liquid Interface (Cambridge 

University Press, London, 1983). 
[13] J. P. Franck and J. Jung, Physica D 23, 259 (1986). 
[14] J. Maurer, P. Bouissou, B. Perrin and P. Tabeling, Eu- 

rophys. Lett. 8, 67 (1989). 
[15] M. Ben Amar and Y. Pomeau, Europhys. Lett. 6, 609 
(1988). 

[16] M. Adda Bedia and V. Hakim, J. Phys. I France 4, 383 

(1994) . 

[17] M. Adda Bedia and M. Ben Amar, Phys. Rev. E 51, 1268 

(1995) . 

[18] D. A. Porter and K. E. Easterling, Transformations in 
Metals and Alloys, (Chapman and Hall, London, 1992), 
p. 114. 

[19] P. W. Voorhees, S. R. Coriell, G. B. McFadden, and R. 

F. Sekerka, J. Cryst. Growth 67, 425 (1984). 
[20] G. P. Ivantsov, Dokl. Akad. Nauk SSSR 58, 567 (1947). 
[21] E. Brener, Phys. Rev. Lett. 71, 3653 (1993). 
[22] A. Karma, Branching in Nature, Les Houches, Vol. 13, 

edited by V. Fleury, J. F. Gouyet and M. Leonetti (EDP 

Scinces/Springer-Verlag, 2001), Chapter XI, p. 365-401. 
[23] C. Zener, J. Appl. Phys. 20, 950 (1949). 
[24] T. Ihle and H. MiiUer-Krumbhaar, Phys. Rev. E 49, 2972 

(1994). 

[25] G. E. Nash and M. E. GUcksman, Acta Metall. 22, 1283 



9 



(1974). 

[26] A. Barbicri and J. S. Langor, Phys. Rev. A 39, 5314 
(1989). 

[27] P. Pelce and Y. Pomeau, Studies in Appl Math 74, 245 
(1986). 

[28] L.M.Fabictti and R.Trivedi, Metallurgical Transactions 
22 A, 1249 (1991); L. M. Fabietti and R. Trivedi, Journal 
of Crystal Growth 162, 185-197 (1997). 

[29] N. Dey and J. A. Sekhar, Acta Metall. Mater. 41, 409 
(1993). 

[30] B. BiUia, D. Boniclli, N. Bcrgeon, Y. Dabo, R. Guerin, 
H. Jamgotchian, H. Nguyen Thi, and P. Voge, Proc. In- 
ternational Conference on Solidification Science and Pro- 
cessing : Outlook for the 21st Century, edited by B.K. 
Dhindaw, B.S. Murty, S. Sen (Science Publishers, En- 
field, NH, USA, 2001) p. 25 



Tables: 

TABLE I. Comparison of the facet length (first column) 
and corner radius (second column) of the equilibrium shape 
for 6 = 1.0 and = 7r/200. Analytical values are compared 
with estimates from phase-field data analyzed with two dif- 
ferent interpolation methods. 



Xr Xl 


(Xt — Xr)V2 










1.950 


1.034 


centered diff. 


approx. [Eq. (40)] 


2.023 


0.983 


one 


-sided approx. [Eqs. (41,42)] 


2.031 


0.978 




analytical 




TABLE II. Convert 


jence of the steady-state tip velocity 


with decreasing ratio Wo /do 


= (Dto) / {aia2Wo) and lattice 


parameter 


h/Wi,). The other parameters are kept constant, 


A = 0.55, S = 1.0, 6»o 


— 7r/200, and the kinetic coefficient is 


set to zero 














Wo /do 


h/Wo 


At/ro 


Vtq/Wo 


Vdo/D 
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5.42 


0.4 


0.008 


0.158 


0.0097 


4 


7.22 






0.285 


0.0098 


5 


9.03 






0.402 


0.0089 


4 


7.22 


0.2 


0.002 


0.285 


0.0098 
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0.008 


0.285 


0.0098 






0.6 


0.018 


0.284 


0.0098 
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0.032 


0.280 


0.0096 






1.0 
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0.275 


0.0094 




e 

FIG. 1. Three functions of the interface orientation an- 
gle 6 represented for a cusp amplitude S = 1.0. Solid line: 
anisotropy function with sharp cusps, f{0). Dashed line: 
smoothed anisotropy function, fs{d). Dotted line: dimen- 
sionless stiffness, fs{0) + fs{0). A large smoothing angle, 
00 = tt/IO, is used here to make the difference between / and 
fs visible. 
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FIG. 2. Analytical equilibrium shape for anisotropy S = 1.0 
and smoothing angle 6o = 7r/200. Pr is the rightmost point 

of the "facet" , Pt the intersection point with the x = y line, 

and Pc the center of the arc circle PrPt- 
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FIG. 3. Comparison of the analytical equilibrium shape 
(line) with the phase-field solid-liquid interface (dots) for 
6 = 0.2, 0.5, 1.0 (from left to right) and Oo = 



x/d 





FIG. 5. Time evolution of a faceted dendrite for undercool- 
ing A = 0.55, smoothing angle do = •7r/200, and anisotropy 
5 = 1.0. The time interval between two succesive curves is 
lOOro. 
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FIG. 7. Variations with the smoothing angle 6o of R/do 
(circles), 0.25(A/do) (squares) and 50{Vto/Wo) (diamonds) 
for anisotropy 5 = 1.0. The continuous lines axe linear fits to 
the data points. 



FIG. 9. Same needle dendrite as in Fig. 5 after rotation, 
translation, and normalization by p = 2Dp/V (solid line). 
Also shown is the corresponding Ivantsov parabola z = 
(dashed line). Note that here, as in Fig. 17, z represents the 
main axis of the dendrite, x the axis perpendicular to z, and 
that both X and z coordinates are normalized by p. 
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FIG. 8. Variations with the cusp amplitude S of R/do 
(circles) and A/do (squares). The continuous lines are 
least-squares fits of the data points to the laws given in 
Eqs. (50, 51). 



FIG. 10. Plots of A/p and R/p for A = 0.55, obtained from 

phase-field sinmlations (squares and circles, respectively) and 
predicted by the sharp-interface theory of section VI (solid 
and dashed lines, respectively). 
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FIG. 11. Tip selection parameter C = Ap^V/Ddo versus 
cusp amplitude 5 for A — 0.55, extracted from the phase-field 
simulations (solid circles) and predicted by the approximate 
sharp- interface analytical theory of section VI (solid- line). 
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FIG. 12. Log-Log plot of the dimensionlcss facet length 
A/do as a function of the dimensionlcss tip velocity Vdo/D 
for a cusp amplitude S — 1.0. The straight line is a linear fit 
to the data points giving a slope —0.54 ± 0.04. 




FIG. 13. Time evolution of a double needle dendrite for 
undercooling A = 0.60 and anisotropy 5 = 1.0 (time interval 
between two contours: lOOro). 
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FIG. 14. Comparison between one branch of a double nee- 
dle dendrite (thick line) and a single needle dendrite (thin 
line). The first curve is shifted to superimpose the two tips. 
Time t = 1500to, undercooling A = 0.60 and anisotropy 
6 = 1.0. 
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FIG. 15. Tip velocity as a function of time for the double 
needle dendrite shown in Fig. 13. The dashed horizontal line 
gives the stationary velocity obtained for the single needle 
dendrite with the same parameters. 
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FIG. 16. Variation of dimensionless tip velocity Vdo/D 
with cusp amplitude 5 obtained from the phase-field simu- 
lations for A = 0.55 (solid circles) and using the equation 
V — 2Dp/p, where p(A) is the Peelet number predicted by 
the Ivantsov relation, and p = {A+R)/y^, where A and R are 
the facet length and tip radius obtained from the phase-field 
shapes (empty circles). 
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FIG. 17. Comparison of the phase-field needle dendrite of 
Fig. 9 (solid line) with the faceted needle crystal shape as- 
sumed in the analytical theory with a circular tip and a 
parabolic tail (dashed line). Both curves are for the same 
parameter value R/{\/2p) = 0.25. 
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